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Abstract 

We analyze the Allan Variance estimator as the combination of Discrete-Time linear filters. We apply this analysis to the 
different variants of the Allan Variance: the Overlapping Allan Variance, the Modified Allan variance, the Hadamard Variance 
and the Overlapping Hadamard variance. Based on this analysis we present a new method to compute a new estimator of the 
Allan Variance and its variants in the frequency domain. We show that the proposed frequency domain equations are equivalent 
to extending the data by periodization in the time domain. Like the Total Variance [1], which is based on extending the data 
manually in the time domain, our frequency domain variances estimators have better statistics than the estimators of the classical 
variances in the time domain. We demonstrate that the previous well-know equation that relates the Allan Variance to the Power 
Spectrum Density (PSD) of continuous-time signals is not valid for real world discrete-time measurements and we propose a new 
equation that relates the Allan Variance to the PSD of the discrete-time signals and that allows to compute the Allan variance 
and its different variants in the frequency domain . 



A. Makdissi is with ELA MEDICAL (SORIN Group), C.A. La Boursidiere, F-92357 Le Plesis Robinson, http://www.alamath.com, 
F. Vernotte is with Institut Utinam, UMR CNRS 6213, Observatoire de Besanfon, Universite de Franche-Comte, 41 bis avenue de I'observatoire, BP 1615, 
F-25010 Besanfon cedex. 

E. De Clercq is with LNE-SYRTE, UMR CNRS 8630, Observatoire de Palis, 61 avenue de I'observatoire, F-75014 Palis 



Stability Variances: A filter Approach. 



I. Introduction 

The Allan Variance [2] and other frequency stability vari- 
ances [3], [4], [5], [1] were introduced in order to allow 
characterization and classification of frequency fluctuations 
[6]. One of the goals of these frequency stability variances was 
to overcome the fact that the true variance is mathematically 
undefined in the case of some power law spectrum [6]. 

The stability properties of oscillators and frequency stan- 
dards can be characterized by two ways: the power spectral 
density (PSD) of the phase (or frequency) fluctuations, i. e. 
the energy distribution in the Fourier frequency spectrum; 
or various variances of the frequency fluctuations averaged 
during a given time interval, it is said in the time domain. The 
power spectral density of frequency fluctuations is of great 
importance because it carries more information than the time 
domain frequency stability variances and provides an unam- 
biguous identification of the noise process encountered in real 
oscillators. PSD are the preferred tool in several applications 
such as telecommunications or frequency synthesis. Stability 
variances are most used in systems in which time measure- 
ments are involved, or for very low Fourier frequencies. Each 
one of these tools corresponds to a specific instrumentation, 
spectrum analyzers for frequency-domain measurements, and 
digital counters for time domain measurements. Although 
there is a separation between measurements methods, use 
and sometimes user's community of these two parameters, 
time-domain and frequency-domain parameters naturally are 
not independent. The true variance for example can be the- 
oretically deduced from the PSD by an integral relationship. 
The true variance cry of a zero-mean continuous-time signal 
Y{t) is defined for stationary signals as the value of the 
autocorrelation function Ry (r) = E[Y (t) Y {t + r)] for 
T ~ (where E is the mathematical expectation operator) 
[7]. This statistical definition of the autocorrelation is related 
to the time-averge of the product Y (t) Y {t + t) if the signal 
is correlation-ergodic [8] by: 
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Ry{T)=\\m^^ I Y{t + r)Y{t)dt 



(1) 



The definition of the two-sided Power Spectral Density 
(PSD) Sy^ (/) of the signal Y is related to Autocorrelation 
function by the Fourier Transform and its inverse by [7]: 



and 



5?^(/)= / i?y (r) e-^2./r^^ 
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(2) 



(3) 



The two-sided PSD is a positive (S'^^if) > 0) and a 
symetric function in / (S'^^{f) = S'^^{—f) ). In frequency 
metrology, the single-sided Power Spectral Density Sy{f) has 
been historically utilized. It is related to the two-sided PSD 
by : 
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For power-law spectrum signals, the PSD is expressed as 
Sy (/) — haf^ [6]. The a integer value may vary from -4 to 
+2 in common clocks frequency fluctuation signals [9]. The 
true variance is defined then as [6]: 



Ry (0) - / Sy if) df 



Krdf. 



(5) 



We can notice easily that for integer a < 0, lim/^o/" 
diverges and then the integral in ^ is infinite. 

The intent of this paper is to explore the relationship be- 
tween stability variances and the PSD using a filter approach. 
This approach allows us to establish new estimators of the 
classical known variances (Allan, Hadamard) in the frequency 
domain instead of the time domain, especially in the case of 
discrete signals, which are the most current in practice. The 
filter approach analysis is developed in Section II in the general 
case of a difference filter of order n. This approach allows 
us to propose general formulae for the stability variance of 
continuous-time signals. The well known frequency stability 
variances like (AVAR, MODAVAR, HADAMARD) are special 
cases of the proposed formula for n=l and n=2. As in 
practical application the signals are not continuous because the 
measurement instruments are read at discrete periodic instants, 
the filter approach is then extended in Section III to discrete- 
time signals. New estimators of the classical variances in the 
frequency domain are proposed which are different from a 
simple discretization of the integral of the continuous-time 
equations. The proposed discrete-time variances are based 
on the fact that filtering in the discrete frequency domain 
is equivalent to a periodization in the time domain. This 
periodization makes our proposed variances estimator have a 
better statistics than the classical estimators. In Section IV 
we present the theoretical calculation of the equivalent degree 
of freedom of the new proposed frequency domain variances 
estimators. Finally, these estimators, the overlapping Allan 
variance (OAVAR), the Hadamard variance (HVAR), and the 
modified Allan variance (MAVAR), are compared in Section 
V to the same estimators in the time-domain using a numerical 
simulation. 



II. CONTINOUS-TIME SIGNALS 

A. Characterization of long term stability by filtering 

Often, it's desirable to characterize the long term stability of 
clocks. Long term behaviour is determined by the components 
of the PSD at low frequencies (/ tends to zero). In order to 
obtain the long term behaviour we average the signal Y{t) 
and we study the variance of the averaged signal. Let Z{t) the 
signal obtained by averaging the signal Y{t) during a time r. 
We can write then: 



Z{t.T) = - [ Y{t)dt. 

T Jt-T 



(6) 



The signal Z {t, t) could be seen as the output of a moving 
average filter M of length r. The moving average filter 
impulse response m {t, r) is defined by: 



. {t, t) = ^Rect, (< - ^) , 



(7) 



where Rect^ {t) is a centered rectangular windows of width 
B: 



Rects (t) = 



1 for - f < t < f 
otherwise. 



(8) 



Thus, in the time domain, Z{t, r) may be defined as: 

Z{t,T)^m{t,T)*Y{t), (9) 

where denotes the convolution product operator. 

The frequency Response M{f) of this moving average filter 
is given by: 



Mif) 



sin (ttt/) 
f — ' 



(10) 



According to hnear filter properties, the PSD Sz (/) of the 
continuous time signal Z {t, t) is: 



Szif) = \Mif)\'SY{f) 



(11) 



From (|5]l (fTOl l and (fTTT l. the variance of the Z {t, r) signal 
is expressed by: 

sin^ {TTTf ) 



-SviDdf. 



(12) 



It's clear that when Sy (/) = ^a/" the variance a\ (r) 
is not defined for power law with a < because the M{f ) 
filter tends to 1 when / tends to zero. In order to make the 
variance cr| (r) defined when a < we need to introduce an 
additional filter D{f) in series with M{f). The input of the 
new D{f) filter is Z{t,T) and let us call its output U {t,T). 
The variance of the U {t, r) signal is expressed when Y{t) has 
power law spectrum by: 

sin^ {ttt/) 



This processing may seem contradictory in the sense that we 
are looking for the long term behaviour (i.e when / — ^ 0) of 
the signal Y{t) and the proposed processing introduces in the 
same time a filter D{f) that eliminates to a certain extent the 
components of the PSD of Y{t) at / = 0. In fact, even if the 
introduced processing may cancel the component of Sy (/) at 
f — and hence makes tend afj (r) to when r — > oo, such a 
processing allows to study the asymptotic behaviour of Sy (/) 
when / approaches zero. We will see in the following of this 
paper that this asymptotic behaviour allows to characterize and 
classify the noise signals. 

In order to realize a filter D{f ) with a frequency response 
D (/) ~ f^, the first idea that comes to mind is to use mul- 
tiple continuous time derivations of the signal Z(t, r). Each 
derivation in the time domain is equivalent to a multiplication 
by in the PSD domain. 

Let (t, r) the rfi^ derivative of Z {t, t) defined by: 



C/(") (t, ' 



d"Z(i,T) 



The PSD S*^"^ (/) of {t, r) is given by: 



2„ sin^ (ttt/) 



and its variance cr^ (t, n) is equal to: 



(14) 



Sy if) , 
(15) 



4 (r, n) = r {2nfr "^^f^Sy if) df. (16) 
Jo (ttt/) 

A perfect continuous derivation in the time domain has 
a linear frequency response for all the frequencies. Such a 
derivation is impossible to realize and is often approximated 
by a filter D that have the same frequency response in the 
vicinity of / = 0. The simplest filter that approximates a 
derivation is the simple time difference filter defined by its 
impulse response: 

dit,T) ^S{t)-S(t-T). (17) 

Its Fourier Transform D{f) is given by: 

£>(/) = !- e-^'^'^f = 2i sin (ttt/) e-'^'^f . (18) 

When cascading n simple difference filters d{t,T), we ob- 
tain a n-order difference filter. Its impulse response d*^"^ {t, r) 
is given by: 



(19) 



fc=0 



Obviously, the variance cr^ (r) becomes defined if 
lim/^o |£) (/)|^ /" is defined. This means that D{f) must 
be of the form when / ^ with (3 > ~a/2. In common 
clock noise with —4 < a < +2 the filter D{f) must verifies 
D (/) cx approximately for sufficiently small / in order to 
make (r) defined for the seven common clock noises. 



D {f)\^ f"df (13) where C„ in the above equation is the binomial coefficient 
defined by C^^ — f,^(^^Lky. '^^^ denotes the factorial of n. 
According to equation ( fTSl ), the frequency response 
(/) of the (t, r) filter is given by: 

£)(") (/) = {D (/))" = 2"r sin" (ttt/) e"*™^^. (20) 

We choose to normalize this filter in such a way that it does 
not modify the variance of a white noise processed by it. The 
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normalization factor c„ is given by the square root of the sum 



of the squares of the coefficients (— 1)'^ C^: 
2 4'T(n + l/2) 

fc=0 



n 



n 
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C^n- (21) 



The output (t, r) of the normalized filter d^") (i, t) /c„ 
is given by: 



Cji L 

The variance of t/^^^^ (t, r) is expressed by: 



(22) 



o'a ('r)(„) = — 



D^^Hf)M{f) Svim (23) 



or, equivalently, using equations ( fTOl i and ( |20] |: 



sm 



2n+2 



(24) 



The convergence domain of this variance is given by a > 
— (n + 2). For positive a values we must introduce a high 
cut-off frequency as the upper limit of the integration in order 
to insure the convergence of cr^ i'^){n)- 

Sometimes it's useful to express the variance afj (''')(„) 
versus the PSD Sx if) of the phase signal X{t) related to 
the frequency fluctuation Y{t) by Y{t) — . Replacing 
Sy if) = [27: ff Sxif) in dH we get: 

22n+2 



sin^"+^(^T/)5x(/)d/. (25) 



Thus, according to the order n of the used difference 
filter (i, r) we obtain different variances with different 
convergence domains (see [4] and [10] for the explicit link 
between n and the convergence). We will see in the next of 
this paper that most of the well-known stability variances are 
special cases of equation (l24l i or dZST l. 

B. The Allan Variance and the Hadamard Variance as filters 

When the order n of the filter d*^"-' (t, t) is equal to one, 
c\ — 2 and, from (l24l l. we obtain the Allan Variance defined 
by: 



sin"^ (ttt/) 



Sy if) df. (26) 



The Allan Variance is noted ct^ (t) in the literature but it's 
the true variance of U'-^^ (t, t), a version of Y(t) processed by 
filters M and D. 

Equation (l26T i shows that the Allan variance is defined for 
power law spectrum with a values from -2 to 0. For a > 0, 
the Allan variance does not converge unless a high cut-off 
frequency fh is taken into account. Moreover the asymptotic 
behaviour of afj (r) is similar for the White Phase noise {a = 
2) and Flicker Phase Noise {a — 1) (see table IJl. For power 
law with a = —3 and a = —4 the Allan variance is undefined 
(unless a low cut-off frequency is taken into account). 

When the order n of the filter d'"-* (i, t) is equal to 2,c\ — & 
and we obtain the three sample Hadamard variance [11] also 



called the Picinbono variance [12]. From (l24l i. this variance is 
defined by: 

sin^ (ttt/) 

(2) = 3 ■ 
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ir) 



-SYif)df. (27) 



This equation shows that the Hadamard variance is defined 
for law power spectrum with integer a values between -4 and 
0. As previously explained, a high cut-off frequency fh is 
necessary for a > in order to ensure convergence of the 
integral dZTl i when / — > 00. 

Table |T] shows the values of Allan variance [6] and 
Hadamard variance [11], [12], [13] for power law spectra. 
The results reported in this table if a > are only valid 
forT» 1/(2^A). 

Because D if) cx /" in the vicinity of zero we can 
say that _D-filtering is equivalent to high-pass filtering. The 
combination of the low pass filter with the high-pass 

filter Dif) forms a band-pass filter G(/). We will see in the 
following of this paper that all the stability variances could 
be expressed as the variance of output of band-pass filters 
applied to the signal under study Yit). When varying t we 
obtain different band-pass filters (a filter bank) with different 
bandwidths. This analysis is similar to the multi-resolution 
wavelet analysis [10] and the special case of the Allan variance 
filter is nothing else but the Haar wavelet basis function [14]. 

It's worth recalling that equation ( |26] | is valid only for con- 
tinuous time signal and filters. This equation gives a theoretical 
definition of the Allan Variance of the continuous signal Yit) 
and cann't be used to compute the Allan variance unless the 
formal expression of the PSD S'y(/) is a known function. In 
real world application signals are collected at discrete instants 
and the above M and D filters are unrealizable for big values 
of T especially when r duration may last for months and years. 
In the next section we analyse the stability variances in the 
case of discrete-time signals. 

III. Discrete-time variances 

In real world applications, measurement instruments are 
read at discrete periodic instants. Let T be the period of the 
reading cycle. We suppose that the instrument measures the 
mean value during this cycle without dead time. We have then 
a discrete time series or signal given by: 

1 f''^ 

yk = T^ Yit)dt. (28) 

J J{k-l)T 

The time-series of a finite length is converted to digital 
numbers and is studied in order to characterize and classify the 
continuous time signal Yit). The PSD Sy if) of the discrete- 
time signal is periodic with a period fs = 1/T and is related 
to the PSD Sy if) of the continuous signal Yit) by: 

Sy if) = TfZ^SYif - nfs) , 7-7^. (29) 



VTif-nf,)X 



We notice from equation ( l29l l that the PSD Syif) is equal 
to SY{f)/T when / ^ because all the terms in the sum are 
null (sin (nvr) — 0) except the term for n — Q. We conclude 
that we can study the long term behaviour of the continuous 
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TABLE I 

Allan and Hadamard variances for power law spectra. 7 ss 0.577216 is the Euler's constant and /,i is the high cut-off frequency 

FOR NOISE WITH O > 0. 
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signal Y{t) by using the discrete time series yk- We can show 
without difficulty that in the presence of a dead-time (sampling 
period larger than the averaging period) we have an aliasing 
phenomenon even for / ^ 0. 

In some applications it's possible to eliminate or reduce 
the aliasing phenomenon by using a low pass filter inside 
the measurement instrument in front of the moving average 
operation. 

For frequencies varying between and /s/2, we can expect 
that the PSD Sy{f) of the discrete sequence yk is nearly equal 
to SyiD/T, at least in the case of a white noise, because 
averaging during a time T and then sampling with a period 
T preserve most of the information contained in the signal 
Y{t), since the averaging can be considered as a non perfect 
anti-aliasing low pass filter. 

For power-law spectrum the sum in equation (|29] | can be 
expressed formally for a < 0. For a > 0, we must introduce 
a high cut-off frequency fh- Table HIl shows the expression of 
Sy{f) for some negative a values when / varies between 
and /s/2. The formulae in Table HIl relating the PSD of the 
sampled signal to the PSD of the continuous signal were never 
published before to our best knowledge. 

A Taylor expansion of Sy{f) when / tends to zero (a < 0) 
gives (see table HHi: 



5. 



y if) — rp 



SY{f) + A[f) 
T 



(30) 

We call A{f) = —ha ^ /"^^ the aliasing term for integer 
a < 0. It depends on the sampling period T and is null for 
white noise (a = 0). At long term, the dominant component 



TABLE II 

The PSD Si, (/) OF THE SAMPLED TIME SERIES J/fc WHEN Y(t) HAS A 
POWER LAW SPECTRUM 5y(/) =haf°'. VC^i 2^) IS THE POLYGAMMA 
FUNCTION DEFINED BY 1p{n,x) = { — l)"+^n\ X^fc^o ^^^^ ^^)"^^ ■ 



TSyjf) 



ho 



[l-T3/3v(2,l + /r)] sm2 (tt/T) 



7r2T2/3 



7r2T2 [2 + cos (27r/T)] 
'~ sin2 (tt/T) 



[12 - r5/5^ (4, 1 + /T)] sin2 (tt/T) 
127r2T2/5 



TT-^T* [33 + 26 cos (27r/T) + cos (47r/T)] 



60 



sin* (tt/T) 



in ( l30b is SY{f)/T and the aUasing is negligible. For short 
term (/ — > ^) the aliasing term varies as T^" and increases 
when the sampling period grows. 

For a > 0, the aliasing term depends also on the high cut- 
off frequency and varies as p whatever the value of a. This 
means that the study of the stability variance of yk for power- 
law spectra with a > 2 does not allow to study the behaviour 
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of Y{t) because the aliasing term is dominant when / tends 
to zero [15], [16]. 

The variance cr^ of the discrete time series yk is related to 
its periodic PSD Sy{f) by [8]: 



T 



' 2T 



(31) 



Equation dSOl l relates the PSD of the measured discrete time 
signal Uk (after averaging without dead-time) to the PSD of 
the continuous-time signal Y{t). Equation ( [3T| ) relates the 
variance to the PSD of the discrete-time signal. Combining 
theses two equations and using an approach similar to that 
presented in paragraph ( III-AI ) in the case of general difference 
filter of order n for the continuous-time signals, allows us 
to define a general stability variance for discrete-time signal 
similar to that of equation (|24]) for continuous-time signals. 

In the case of a frequency fluctuation sequence, the time 
series yk could be related to the time error samples X{t) by: 



T 



kT 



(fe-l)T 



dX{t) 
dt 



dt 



X [kT] -X[{k^ 1) T] 



(32) 



Sometimes, it's difficult to realize experimentally the mea- 
surement of yk according to equation ( l28T l by averaging and 
recording y^ without dead-time. If the time error data X{t) are 
measurable it is always possible to sample them and compute 
yk according to equation ( |32] | without dead-time. 

In order to simplify notations, we suppose, without loss in 
generality, that T is equal to 1 in the following of the paper 
Then, integration in equation ( |3T1 ) is done over the interval 
[—1/2, 1/2] and equation ( |32] | could be written, by denoting 
Xk^X [kT), as: 



Vk = Xk - Xk^i- 



(33) 



In other terms, the time error sequence Xk could be ob- 
tained from the averaged frequency signal yk by numerical 
integration with a starting point xq ~ 0: 



Xk+i = Xk+ yk- 



(34) 



In order to estimate the (T^(r)(„) variance from the ob- 
served discrete-time series yk we try to realize a discrete 
version Uk of the continuous signal t/'"^ (t, t) defined by 
equation (l22l i by using digital filters similar to the analog 
filters m(t, r) and (i^"^(t,r). Once we have a discrete version 
of U^'^\t,T), we can estimate its variance by computing the 
sample variance of the discrete-time series Uk- 

Following the filter approach used for continuous time 
signals we introduce digital filters in such a way that their 
discrete-time outputs are similar, as much as possible, to 
analog signals in the previous section. 

The moving average filter m{t, r) of length t becomes in 
the discrete domain a rectangular windows of length m — 
t/T. The output Zk of this filter is given by: 



Zk 



^ m — 1 

— Y\ vk- 



(35) 



Obviously, by using (|35] | and (|28] |. we can write: 



Zk 



1 

mT 



kT 



Y (t) dt. (36) 

'(fc-m)T 

Equation (|36] | shows that averaging m values of the signal 
yk is equivalent to using an instrument with an averaging time 
T = mT. This may let us think wrongly that the PSD Sz (/), 
of the discrete time series Zk could be obtained directly from 
equation STT[ by replacing t = mT. 

In fact, Zk being discrete, its PSD is periodic and contains 
aliasing terms. The PSD Sz (/) of the discrete time series Zk 
is related to the PSD Sy (/) of the continuous signal Y(t) by: 



f 2^ Sy [f - nfs) — — zrrz ttt^-- (37) 



T 



[nmT (/ - nfs 



In order to relate the The PSD Sz{f) of the averaged 
discrete time series Zk to the PSD of the sampled signal 
yk we compute the Fourier Transform M* [T) of the digital 
filter TTifc where is a normalized frequency for the discrete 
time signals: T — f ■ T . The impulse response of this filter 
is TOfc — —TTjn (k), where tt^ (k), is a discrete rectangular 
window of length m with all its coefficients equal to 1 . This 
impulse response is obtained from m {t, t) by sampling it with 
a sampling period T. The Fourier Transform M* {T) is then: 



m — 1 



M*{T) = 



m ^ — ^ 



fe=0 



-2i7^kj^ ^ 1 1 - exp (-2ii7rJ^m) 
ml — exp (— 2j7rJ-") 



1 sin (ttttiJ-") 
m sin (ttJF) 



exp [— zttJF (to — 1) 



(38) 



We can notice that the frequency response of the discrete 
moving average filter of equation dSSl) is different from that 
of the continuous moving average filter of equation dTOl ) when 
replacing t by mT. 

As for the continuous time signals, this M filter is not 
sufficient to ensure the convergence of the variance for power 
law spectrum signals with a < 0. Therefore, we introduce 
a digital version of the continuous D filter by choosing an 
impulse response dk as: 



dk — (Sk — ^fc-m) 



(39) 



where Sk is a Dirac impulse of unity amplitude. 

As for the discrete time filter mk, the discrete filter dk is 
obtained by sampling d {t, r) of equation ([TTT i. 

The frequency response D* {T) of the digital filter dk is 
identical to that of the continuous filter D{!F): 



D* {T) = l-e 



— 2inJ-'m 



2ism{7rTm)e-'''^"'. (40) 

(n) 

When using n difference filters we get the digital filter rf^, 
by sampling the continuous time filter d*^") (t, r) of equation 
03: 



p=0 



(41) 



n=0 



This impulse response could be obtained also by a digital 
convolution (denoted by (g) in the following) of the filter dk 
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in equation ( l39b with itself n times. The frequency response 
£)*(«) (jr) Qf tijg fljtgj (^C'l) according to (gOll, given by: 



D<'^\T) = [L>* (J^)]" = 2"i" sin" (ttJ^to) e" 



■imrj-'m 



(42) 



If we use the same normalization factor c„ as the ones of 
equation (1211 1. the output Uk of the normalized filter d^f^ jcn 
is given by: 



(43) 



According to equations ( |43] |. dlTT i. dSST l and (l42l l. the true 
variance cr^ (m) of the discrete signal is related to the 
PSD 5*1, (/) of the discrete signal by: 



(™) 



ct.ra'- 



1/2 gi^2n+2 



{irfm) 



2 2 1 • 2, n (44) 



1/2 



Equation (l44l i defines a stability true variance of discrete- 
time signals in the general case. Percival proposed in [17] an 
identical formula to that obtained in (l44l l when n = 1 in the 
case of Allan variance. 

Comparing this expression to equation (|24] | we can notice 
that the denominator in (l44l l is sin (vr/) while that of 
equation (l24l) is (ttt/)^. We have shown in equation (l30l l that 
Sy{f) « S'y(/)/r. This difference bewteen equations (l24l l 
and (l44l l may let us think that the true variance cr„ (m) of 
Mfc is different from the variance afj (''')(„■) of the continuous 
signal f/^"^ (iiT)- Appendix I VII show a mathematical demon- 
stration of the equivalence of the discrete-time variance and 
the continuous-time variance. 

The above discret-time variance can be written versus the 
PSD of the discrete-time error samples Xk- Using equation 
(l35T l and ( |33] | we can write: 



mzk = Xk - Xk- 



^k 



(45) 



Using this expression in ( l43T l we can express Uk in terms 
of the phase measurement Xk under the simple form: 



Uk = 



2 C„ V 



'E)Xk 



(46) 



It's clear that equation (|46] | is simpler than equation (|43] | in 
terms of computation complexity because the filter d!]f'®mk 
of equation ( |43] ) must be computed for each m value while 
the coefficients of the filter rfjT'^^^ of equation ( |46] | do not 
depend on the averaging factor m. 

According to equations (|46] |, (ISTT i and (l42T i. the true variance 
cr^ (m) of the discrete signal Uk is related to the PSD Sx (/) 
of the discrete signal Xk by: 



22n+2 



-1/2 



-1/2 



,2n+2 



(TT/m) 5, (47) 



This equation shows that the transition from the stability 
variance of the continuous-time signal X{t) given by equation 
(I25I 1 to the stability variance of discrete-time signal Xk is done 
very simply. 



A. Estimation of the Stability Variances of the Discrete-Time 
Signals 

In order to estimate the variances presented in the last 
section we use the sample variance of the zero mean discrete 
signal Uk- 

1 



N-l 

N ^ 

k=0 



Uk 



(48) 



where N is the length of the time series Uk- 

When the signal Uk is obtained by filtering a signal of length 
L using a filter of length p, we must consider in (l48T l only 
N — L ~ p + 1 unambiguous samples of Uk- 

Let Uk be the Discrete Fourier Transform (DFT) of the 
discrete signal Uk defined by: 



Un 



N-l 

E 

fc=0 



UfcC 



, ne{0,---,7V_l}. (49) 



The sample variance can be related to the DFT series using 
the discrete Parseval's theorem: 



E 

fc=0 



N-l 



Uk\ 



^Ei^^i^ 



(50) 



fe=0 



The Uk coefficients for N/2 < k < N represent the 
negative frequencies. In the case of a real signal Uk, the 
coefficients Uk are symmetrical around P = [-^^^j- We 
define a "one-sided" set of DFT coefficients Uk by: 

Hi 

Uk = Uk, 0<fc<P-l 

( Up if N is odd 
Up^ { Up_ 
V2 

The Parseval's theorem could be written then: 

p 



Uk 



(51) 



if N is even. 



N-l 

E Kl' 

fc=0 



fe=0 



(52) 



According to equations ( l43T l and (l45b . the DFT coefficients 
Uk of the time series Uk are related to that of Xk and yk 
by: 



Uk 



N 



^,Yk = f^) Xk. 

N I m Cn \N J 

(53) 

The transition from equation (03]) to the first part of the 
above equation is valid under the assumption that discrete- 
time signals are N-periodic. This means that the sample 
variance in the frequency domain is equivalent to the sample 
variance in the time-domain applied to an extended version 
(by periodization) of the discrete-time signal. The first part of 
above equaUty gives when using (l38T l. (l42T i. ( |52] | and ( l48T i: 



22n+l 



E- 

k=0 



-,2ri+2 



N J 



Yk 



(54) 



To our knowledge, this is the first time that a relation be- 
tween the sample variance estimator of the frequency stability 
and the DFT of discrete time series yk is established. It's 
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worth recalling that this equation is not a direct approximation 
to compute the generic variance expression of equation (l24l l 
by discretization in the frequency domain as was proposed in 
[18] but it is the variance, according to the Parseval's theorem 
( |52] |. of a signal yk filtered in the frequency domain . 

Some works [19] have shown that using the numerical 
integration in (l24l l to estimate the Allan variance (n = 1) leads 
to a biased estimator regarding the classical Allan variance 
sample estimator We will show at the end of this paper that 
our formula ( |54] | gives results which are nearly identical to 
the classical sample estimators. 

In fact, if we can consider that Y{t) is band-limited to 
/max — Yr ^^^^ ^® approximate the integral in equation 
(l24l l in the Riemann sense by replacing the integration by the 
sum of the surfaces of rectangles of width at discrete 
frequencies fk = 



p 

y 

" fc=0 



sm' 



2n+2 fTTkrn 



fc2 



\NT J ' 
(55) 

where Syif) is an estimator of the PSD Syif)- If we use 
2T\Yk\'^/N as an estimator of Sy [fk) then equation 
becomes : 



D"'(l) 




M(m) 




Z)l"(m) 







M(m) 



Fig. 1. The processing chain of the stability variances, M(m) is a moving 
average filter of length m. D^"\m) is a difference filter of order n and lag 
m. t m is the decimation by a factor m operator. 



interesting result could be written as: 



(60) 



In other words, the sample variance of equation (|54] | is equal 
to the integral of equation ( l24l i for a band-limited Y(t) when 
evaluated in the Riemann sense over the interval / e [0, 1/2T] 
by using the periodogram of yt as an estimator of the PSD 
Sy (/) of Y(t) according to equation (|59] |. 

The second part of equation (|53] | gives when using (l42l . 
^ and (|48]i: 



22n+l 



P 



sm 



2ri+2 



\ N ) 



\Yk\' 



(56) 



fc=0 



It's clear that equations (l56T l and ( l54l i are different. This 
difference could by explained by the fact that equation (|24] | 
is given versus Sy (/) which is not observable directly 
while equation i55[ use lYfep, an estimator of the PSD of 
the averaged and sampled version of Y{t). In other words, 
averaging according to equation ( |32| | is considered when using 
iFfeP in equation dSUl while Sy (/) in equations (OHi and (l55T l 
is considered before averaging according to equation (|6]l. 

In order to relate equation (|55] | to equation (l24l i we suppose 
that 5y (/) is band limited. In this case, there is no aliasing 
in equation ( |29] l and it could be written: 



for < / < — . (57) 



1 sin^ (nTf) 

Sy if) - ^Sy if) 



The DFT coefficients Yk could be considered as an estimator 
of the PSD Sy (/) of the discrete signal yk at discrete 
frequencies /&: 



Sv{fk) = ^\Yk\ 



for < fc < P. 



(58) 



This equation is known in the literature as the periodogram 
spectrum estimator The factor 2 in ( |58] | is due to the fact that 
the PSD Sy{f) is one-sided. 

Replacing equations ( |58] | in ( fSTl l. we get an estimator of the 
PSD Sy{f) of the band-limited continuous time signal Y{t): 



Sy 



k 



2T{TTky 



m 



Yk 



(59) 



Using this expression in equation (1551 1 leads to an expres- 
sion identical to the sample variance of equation (l54l i. This 



l2ri+3 



E 

fc=0 



sm 



2n+2 



/ Trkm\ 








\ N J 





(61) 



When X(t) is band-limited, equation ( |6TI ) can be obtained 
directly from equation dZSb using a Riemann sum and replac- 
ing Sx if) by the periodogram of the discrete signal Xk- 

In the following of this paper we express the different 
stability variances in the discrete time using the signal Uk- 
Figure [U shows the different filters involved in the computation 
of theses stability variances. 

B. The Overlapping Allan Variance (OAVAR) 

This is a special case of the above processing when the order 
of the difference filter n is equal to one. The normalization 
factor ci is given by equation (l2Tl i and is equal to \/2. The 
filter d^^^ of equation (ffSI l is equal to 5k — 25k+m + <5fc+2m- 
The signal Uk is given by: 

1 



Uk 



ixk+2m — 2Xfc+ 



(62) 



(2) 



Let N the length of the discrete time series yk- The u^, 
filter length is 2m and the output Uk length is — 2m + 1. 
According to equation ( |48] ), the sample variance of Uk is: 



OAVAR(m) = al (m) (63) 

^ N-2m 
7: — T y ixk+2ni - 2xk+m + Xk)"^ 

— 2m + 1 ) '^-^ 



2m2 {N - 



k=0 



which is the classical estimator of the Overlapping estimator 
of Allan Variance [20]. 

The computation in (|63] | from Xk requires four additions 
and one multiplication for each term inside the sum. The sum 
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over k requires TV — 2m + 1 addition. The whole computation 
requires roughly 5 x N operation and is linear in N. 

When the available measurement are frequency fluctuations 
( \yk\), it's more efficient (in number of floating point opera- 
tions but not in memory use) to compute the phase signal Xk 
using (|34] | and then use ( l63T l to compute the OAVAR variance 
than to compute Zk from yk and then Uk- 

Replacing n by 1 in equation (|54| | we get an expression 
of the Overlapped Allan Variance versus the one-sided set of 
DFT coefficient Yk of the measurement time series yk by: 

4 ^ sin"' ( 2^Hi") 
OAVAR(m)^ = (m) = ^1 " IT^I^^I'- 

(64) 

The DFT computation complexity is N log(N) when using 
a Fast Fourier Transform (FFT) algorithm. But the most CPU 
consuming in ( l64l l is the computation of the sine trigonometric 
functions inside the sum symbol. It's trivial that the computa- 
tion using equation ( l63T l is more efficient than using equation 
(EUl. 

It's worth recalling that the discrete time formula ( |63] | use 
N — 2m + 1 terms. The largest acceptable m value is iV/2. 
In this case the variance is estimated from one sample only. 
The DFT formula ( |64] | use P terms whatever the m value. 
When m = N jl the half of the sine terms in (|64] | is null. 
The computation of the confidence levels when using equation 
(|64] | has shown that the confidence levels are better than 
that of the discrete time formula of equation (l63T l because 
filtering in frequency domain use all the available samples 
while filtering in the time domain use N minus the filter 
length samples. In fact. Filtering in the DFT domain is done 
by multiplication of the DFT. This multiplication is equivalent 
to circular convolution in the time domain. Circular or cyclic 
convolution of two signal of length N is equivalent to classical 
sum convolution with indices modulo N . This means that DFT 
formula ( |64| | is equivalent to a kind of Total Variance [1] where 
the series yk is extended by periodic (circular) repetitions. The 
Total Hadamard Variance [9] uses an extended version of yk 
where the extension use a reflected copy of yk- 

C. The "Non Overlapping" Allan Variance (AVAR) 

The "Non overlapping" Allan variance is a special case of 
the classical Allan Variance that doesn't use overlapped values 
when computing ^ u^, in the sample variance of Uk- This 
means that only {N — 2m + \ )/m values are considered when 
forming the sum. 

In other words, the Allan Variance AVAR is obtained from 
Uk by a decimation operation of order m (See Figure [TJ. If 
we start the decimation at fc = we can use N/m— 1 values. 
The decimated signal is given by: 

N 

fk = Ukm , < fc < 2. (65) 

m 

Replacing (l65T l in (l63T l we get the non overlapping Allan 
variance as the sample variance of r^: 

AVAR(m) = (72(m) (66) 

~ 2m{N-m) 2^k=0 \X(k+2)7yi ^X{k+l)7yi + Xkm) 



It's obvious that the AVAR requires less computation than 
the OAVAR. In fact, for each m value there are N/m — 1 
terms. The largest acceptable m value is Nl2. In this case 
the sample variance is estimated from one sample only. The 
confidence levels for AVAR and OAVAR are equals for m — 1 
and m = N/2. Values of m between m = 1 and m — N/2 
give a better confidence levels in the OAVAR than in the AVAR 
variance. 

Because the OAVAR confidence levels are globally better 
than those of AVAR , the only interest to use the AVAR instead 
of OAVAR is its computation efficiency. 

Though decimation operation of equation ( l65T l is very sim- 
ple in the time domain it has no interest in the frequency 
domain. In fact, the computation of the DFT coefficient Rk 
of rk versus the DFT coefficients of Uk is given by: 

^ m — 1 

Rn = — U(^„^pN/m) Modulo N. (67) 

fc=0 

Computing the sample variance of rk according to (l52T l 
require an additional loop to compute Rk- For this reason we 
don't propose a formula to compute the AVAR in the frequency 
domain as we did for OAVAR in equation $64\ . 

D. The Modified Allan Variance (MAVAR) 

The modified Allan Variance was introduced [5] to over- 
come the relatively poor discrimination capability of the Allan 
variance against white and flicker phase noise. 

Let 7fc be the signal obtained from Uk by a moving average 
filter M (m) of length m (See Figure [ij: 

7fc = — {uk + Uk+i H h Uk+m-i) ■ (68) 

m 

Using Uk expression from (|62| | in ( l68T l we get: 

V2m'^jk = Xk-\ ^Xk+„i-i 

— 2 (Xfe+„i + • • • + Xk+2m-l) 
+Xk+27n-\ \-Xk+3m-l- (69) 

Using this expression directly to compute jk requires a 
summation loop with 3 * (m + 1) floating point operation. The 
biggest acceptable m value in this equation is m = iV/3. This 
yields a computation complexity of iV^. 

In order to reduce the computation complexity we propose 
a recursive formula. Expressing Ak+i = \/2rr?^k+\ using 
equation ( |69] l we can write: 

^fc+l = Ak + Xfc+Sm - 3a;fc+2m + SXfe+m - Xk (70) 

with a starting value computed using (l69l l with fc = 0. 

Allan [21] already proposed a recursive method in order 
to reduce the computation complexity of the Modified Allan 
Variance without giving the details of the recursive equation. 

The computation complexity of Ak according to ( iTOl i is 
linear in N . 

The length of the time series Uk v& N ~ 2m + 1 and the 
length of the filter M{m) is m. we conclude that the length 

of 7fc is — 3m + 2. 
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The Modified Allan Variance MAVAR is the sample vari- 
ance of 7^: 

-1 



MAVAR(m) = a'i (m) 



1 



2m4 [N -2,7X1 + 2) 



Ar-3m- 

E 

fc=0 



(71) 

The PSD (/) of the discrete time signal 7^ is related 
to that of Uk by: 

\M* if)D* if )M* if)\'Syif) 
2 sin^ (ttto/) 



-Sy if) . 



m"* sin** (tt/) ^ 
The DFT coefficients r„ of the series 7^ are given by 



(72) 



1;, 



N ) 



(73) 

According to the Parseval's equation (l52T i for the series 7^ 
we can express the MAVAR versus the one-sided set of DFT 
coefficients of the measured signal by: 

p 



MAVAR(to)f = ("i) 



\ N ) 



m^N^T^^ sm- 

/c— 



(74) 

As for the OAVAR formula in the frequency domain this 
equation is an estimator of the MAVAR in the frequency 
domain. The main difference with the discrete time formula 
dTTT i is the number of terms involved in the sum: P in the case 
of equation (l74b and N — 2m + 1 in equation dTTT l. 



E. The Overlapping Hadamard Variance (OHVAR) 

This is a special case of the above processing when 
the order of the difference filter n is equal to two. The 
normalization factor 0-2 is given by equation (I2TI) and is 
equal to The filter rf^^' of equation ( l46b is equal to 

Sk - iSk+m + 34+2m - 4+4m- Wc dcnotc gfc = Mfc whcrc 
Uk is given by (|46] | with n = 2: 



1 



9k 



(75) 



(3) 



Let N the length of the discrete time series i/k- The d)^ 
filter length is 3m and the output gk length is — 3m + 1. 

The Overlapping Hadamard Variance is the sample variance 
of gk- 



As for formulas ( |64|) and (f74b . equation dTTl) is a new 
formula that allows to compute the OHAVAR in the frequency 
domain. 

It's clear from equation dTSl l that the Hadamard variance 
estimator in the time domain cancels linear drifts. In fact, if 
Uk — k then Xk = k {k — 1) /2 according to equation ( l34b . 
Replacing this value in equation ( fTSI l leads to gk — whatever 
the value of m. 

F. The Hadamard Variance (HVAR) 

The Hadamard variance is a special case of the Overlapping 
Hadamard Variance that doesn't use overlapped values when 
computing J^dk '^^e sample variance of gk- This means 
that only {N — 3m+l)/m values are considered when forming 
the sum. 

In other words, the Hadamard Variance HVAR is obtained 
from gk by a decimation operation of order m (See Figure 
[T]|. If we start the decimation at fc = we can use N/m — 2 
values. The decimated signal hk is given by: 



hk — Qkn 



N 

0<k< 3. 

m 



(78) 



Replacing ( fTSI l in dTSI ) we get the non overlapping 
Hadamard variance as the sample variance of hk'- 



HVAR(m) 



al (m) 



1 



N/m-3 

E ' 

k=0 



6m (N — 2m) 
-3x(^k+2)m + 3a;(fc+i) 

m 



' - (79) 



As for the Non Overlapping Allan Variance AVAR we don't 
propose a formula in the frequency domain for HVAR because 
the decimation operation doesn't simplify computation in the 
frequency domain as it does in the time-domain. 

IV. Frequency variances Equivalent Degree of 
Freedom 

We can express the frequency-domain variance estimator by 
the general form : 



p 

E 

k=Q 



Hk{n, m) 



Y 



N 



(80) 



OHVAR(m) 



(to) 



Where n is the difference filter order, m is the averaging 
(76)factor and Hk{n,m) is given by : 



1 



1 



3x 



fe+2m 



■ 3a;/c4 



Xk) 



6to^ (N — 3m 

7V-3m 

X ^ (Xk+S 
fe=0 

Replacing n by 2 in equation ( l54l ) we get an expression of 
the Overlapping Hadamard Variance versus the one-sided set 
of DFT coefficient Yk of the measurement time series yk by: 



Hk{n,m) 



■)2n+l 



sm 



2n+2 



\ N ) 



clmm sin2 {^) 
for the non-modified variances and : 

2n+4 / Trkm ^ 



~,2n+l 



Hk{n,m) 



sm 



/ Trkm \ 
\ N J 



clm^N 



sm 



(f) 



(81) 



(82) 



OHVARf (m) = g (m) 



16 



3to2^2 



P - 6 

sm" 



nnik ^ 

N , 



for the modified parlances. 



\Yk\\ 



(77) 



The quantity 



Y 



at discrete frequency values fk 
written as : 



jN is the periodogram P(/) evaluated 
Equation ( [8O] ) can be 



k_ 
' N- 
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(83) 



fe=0 



The periodogram P{f) is an estimator of the PSD Sy{f) 

We estimate the Equivalent Degree of Freedom (edf) of ^l* 
by: 



2 (£;(*))= 



Var {^) 

The mean value E (5*) is given by : 

p 

E{^) = Y,Hk{n,m)E{P{h)) 

fc=0 



(84) 



(85) 



It is well know that the periodogram is a biased estimator 
of the PSD Sy{f ) and that : 

E{PU))^WB{f)®S{f) (86) 
Where Wb (/) is the Bartlett window defined by : 

TV sm (tt/J 

and (g) denotes the circular convolution defined by : 

.1/2 

WB{f)®S{f)= WBi9)S{f~9)d9 (88) 

J -1/2 

It's clear that the periodogram is asymptotically unbiased 
since as N becomes very large Wb (/) approaches an impulse 
in the frequency domain. Then we can write for large : 



E{^) = J2Hk{n,m)S{fk) 



k=0 

and for power law spectrum : 

p 



k 

2N 



(89) 



(90) 



The variance Var (^f) is given by : 



Var{-^) =E{^^) 
E E Hk{n,m)H,{n,m)CoviP{fk),P{M 

k=Oj=0 



(91) 



The covariance of the periodogram is given by : 

Cov{P{h),P{h))= (92) 

Replacing /i hy fk ^ j^- and /2 by fj = ^ in equation 
we get : 



Cov {P{fk),P{f,)) = 

Sy ifk) Sy (/,) (^^^|^j£|^£L_^ 



(93) 



TABLE III 

Computation time in ms of the different stability variances, 
N = 400000 





AVAR 


OAVAR 


MAVAR 


HVAR 


OHVAR 


Time Domain 


16 


47 


78 


16 


47 


Frequency Domain 




265 


265 




265 



TABLE IV 

Computation time in ms of the different stability variances, 
N = 2000000 





AVAR 


OAVAR 


MAVAR 


HVAR 


OHVAR 


Time Domain 


63 


265 


484 


63 


360 


Frequency Domain 




1453 


1500 




1485 



Therefore, the covariance ( l93T l is is seen to go to zero when 
k j . The variance is therefore : 



The edf is, according to ( l84l ). given by : 



p 

Y,Hl{n,m)Sl{h) 

k=0 



edf^ 



2{j:k=oHk{n,m)Sy{h)f 
Ek=oHl{n,m)Smk) 



For power law spectrum we get : 



edf 



(94) 



(95) 



(96) 



With Hk {n, m) given by ( l82b for the modified variances 
and dSTT ) for the non-modified variances. 

V. Time Domain versus Frequency Domain: 
Numerical Results 

We have simulated time series data i/k of length N = 
400000, 2000000 and 65536 for the different power law spec- 
tra for —4 < a < 2. Table |TII] and |IV] show the computation 
time on a personnal computer (pentium IV or equivalent @ 
2.8 GHz) in ms of the different stability variances mentioned 
in this paper The computation time of the FFT was included 
in the computation time of the frequency variances. 

For the computation in the frequency domain we used the 
FFT algorithm of Cooley and Tuckey[22]. The FFT com- 
putation time is 45 ms for N — 400000 and 250 ms for 
N = 2000000. 

We presented in equation ( |54] | a new way to compute the 
different stability variances using the DFT of the data. We 
demonstrated that this equation is equivalent to the equations 
in the time domain with a slight difference in the number of 
samples when computing the sample variance. For example, 
equation (l63l l in the time domain uses only unambiguous 
samples in the sense that a filter of length 2m will produce 
N — 2m + 1 unambiguous output samples when applied to an 
input data of length N. 

In the following we present numerical results of the differ- 
ent frequency domain variances estimators presented in this 
paper. The error bars on the plots were computed using one 
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F-OAVAR edf 




10° 10' 10^ 10° 10^ 

Averaging Factor 



Fig. 2. F-OAVAR edf for three noise types for sequences of length N = 
65536. WHFM for White frequency noise, FLFM for Flicker frequency noise 
and RWFM for Random Walk frequency noise. The continous lines (denoted 
"TH" on the Figure legend) represent the theoretical edf computed by equation 
496) . The symbols (denoted "MC" on the Figure legend) represent the edf 
obtained by Monte Carlo simulation with 1000 trials. 
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Fig. 3. OAVAR computed in the time domain and in the frequency domain 
for a White Frequency Noise sequence of length = 65536. The spectral 
OAVAR estimates were slightly shifted in order to be distinguished from the 
time OAVAR estimates. The dashed continuous line represents the theoretical 
response /io/(2t). 

sigma Chi-squared distribution with an equivalent degree of 
freedom {edf) estimated by making Monte Carlo simulations 
of 1000 trials. 

A. OAVAR 

Figure (|2]) depicts the edf of the Overlapping Allan Variance 
computed in the frequency domain (F-OAVAR) for three noise 
types: a White frequency noise (WHFM), a Flicker frequency 
noise (FLFM) and a Random Walk frequency noise (RWFM). 
It shows a very good agreement between the theoretical edf 
formula of equation ( |96] l and the edf obtained by Monte Carlo 
simulations. 

Figure (O compares the Overlapping Allan variance of a 
white frequency noise sequence computed in the time domain 
and in the frequency domain from relationship (|64] |. No bias 



TABLE V 

Comparison of the equivalent degrees of freedom {edf) of the 

TIME T-OAVAR ESTIMATES, THE SPECTRAL F-OAVAR ESTIMATES AND 

THE Total variance estimates for a White Frequency Noise 

SEQUENCE OF LENGTH A'' = 65536. 



T 


T-OAVAR 


F-OAVAR 


TotVar 


1 


46591 


42297 


45368 


2 


40640 


37232 


34379 


4 


24186 


23639 


22460 


8 


11870 


12338 


11451 


16 


5865 


6786 


6375 


32 


2937 


3255 


2945 


64 


1493 


1515 


1555 


128 


746 


740 


832 


ZJO 


3o3 


372 


414 


512 


199 


194 


215 


1024 


93 


89 


104 


2048 


43 


43 


53 


4096 


20 


22 


26 


8192 


10 


12 


12 


16384 


4 


6.4 


6.2 


32768 


1.0 


3.0 


2.9 



is visible between these computations and the theoretical 
response (less than 1 %). On the other side, the error bars 
of OAVAR computed in the frequency domain are clearly 
smaller as the ones of OAVAR computed in the time domain, 
as expected in section Illl-BI Table [V] shows the equivalent 
degrees of freedom {edf) of the Total Variance and the OAVAR 
estimates in the time domain (T-OAVAR) and in the frequency 
domain (F-OAVAR), assuming a Chi-square statistics [23]. For 
the highest r value (t — N/2), the edf of the spectral estimate 
is 3 times higher than the edf of the time estimate, i.e. the 
spectral estimate is \/3 times more accurate than the time 
estimate. 

Such an advantage is particularly useful for detecting and 
measuring the level of the low frequency noises (e.g. random 
walk FM) sooner as with time variances, i.e. for shorter 
duration. Considering that the edf decreases approximately as 
r^^, an estimator with an edf 3 times higher than another 
one provides a noise level estimation \/3 times sooner than 
the other one (e.g. 7 month instead of 1 year) with the same 
accuracy. 

Figure (|4]i presents a comparaison between the Overlapping 
Allan variance computed in the frequency domain (F-OAVAR) 
and the Total variance for three noise types : WHFM, FLFM 
and RWFM. The upper plot depicts the edf ratio computed 
using Monte Carlo simulations with 1000 trials, we notice 
that the edf of the F-OAVAR and the Total variance are nearly 
identical. The lower plot depicts the bias defined by Bias = 
100 X (1 - ^F-OAVAR/Totvar). The bias of the F-OAVAR 
with respect to the Total variance is less than 10%. 

In the same way, figure (|5]l presents a comparaison between 
the Overlapping Allan variance computed in the frequency 
domain and the classical Overlapping Allan variance computed 
in the time domain. The upper plot shows that the F-OAVAR 
edf is two to three times higher than the edf of the T-OAVAR 
for the higher r value (r = N/2) .The lower plot depicts the 
bias defined by Bias = 100 x (1 - ^F-OAVAR /T-OAVAR) . 

Figure ^ shows the Total Variance, the Overalpping Allan 
Variance computed in the time domain (T-OAVAR) and in the 
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-WHFM 

-FLFM 

-RWFM 




Fig. 4. Comparaison of the F-OAVAR and the Total variance for three noise 
types. The upper plot depicts the edf ratio and the lower plot depicts the bias. 
N = 56536. Results were obtained using Monte Carlo with 1000 trials. 
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Fig. 6. The T-OAVAR, the Total variance and the F-OAVAR for a White 
frequence noise (a = 0) and a Ficker frequency noise (o = — 1) . A linear 
frequency drift was added to the noise sequences of length N = 56536 
(Monte Carlo trials = 1000). 



N 




Fig. 5. Comparaison of the F-OAVAR and the T-OAVAR for three noise pjg 7 Random Walk Frequency Noise sequence: rough (A), drift removed 
types. The upper plot depicts the edf ratio and the lower plot depicts the bias. ^g) ^^d circularized (C) 



N = 56536. Results were obtained using Monte Carlo with 1000 trials. 



frequency domain (F-OAVAR) for a White frequence noise and 
a flicker noise with a linear frequency drift. The added linear 
drifts is equal to D{t) ^ 15t . Like the the Total variance and 
the classical Allan variance, the F-OAVAR does not cancel the 
linear drift. We can notice also that the F-OAVAR for a linear 
drift varies as r , while the Total variance and the T-OAVAR 
vary as t^. 

Unfortunately, the last result shows that the computation of 
OAVAR in the frequency domain presents a severe drawback: 
it is unable to discriminate between a linear frequency drift and 
a /^^ frequency noise (random walk FM). This effect is due 
to the assumption of periodicity of the sequence implicitely 
induced by the use of the FFT algorithm. Figure |7]-A shows 
that connecting the last sample to the first one may induce a 
high edge, altering the variance measurements. So we decided 
to process the frequency deviation sequence with 2 different 
ways: 



• by removing the linear drift of this sequence (see figure 
[Tj-B; let us notice that there is still an edge at the end of 
the sequence). The removed line is estimated by a least 
squares fit of the data sequence to a line. 

• by circularizing the sequence (see figure |7]-C), i.e. by 
removing the linear drift in such a way that the last 
sample of the residuals is equal to the first one. Denoting 
hy D{t) = a ■ t + b the drift we have to substract from 
the sequence, the linear coefficient a is then: 

a = yji-^ (97) 

tN — h 

and the constant term b may be choosen equal to since 
OAVAR is not sensitive to additive constants. 
It is worth recalling that Figure (|7]i shows the side effect 
of periodization (induced by multiplication in the discret 
frequency domain) of a sequence without processing, after a 
line removal, and after circularization. But when computing 
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TABLE VI 

Comparison of the equivalent degrees of freedom of the time 
oavar estimates and the spectral oavar estimates rough, 
after removing a linear drift and after circularizing the 
sequence for a random walk frequency noise sequence of 

LENGTH Af = 65536. 



T 


Time OAVAR 




Spectral OAVAR 






rough 


without drift 


circularized 


1 


68540 


65660 


39 


56735 


2 


35289 


33269 


39 


27589 


4 


15498 


15410 


39 


13009 


8 


7324 


7725 


38 


6392 


16 


3621 


3997 


38 


3258 


32 


1812 


2091 


37 


1737 


64 


900 


1040 


36 


860 


128 


455 


477 


35 


436 


256 


225 


219 


34 


224 


512 


no 


106 


31 


109 


1024 


52 


54 


25 


50 


2048 


25 


28 


17 


23 


4096 


12 


14 


10 


11 


8192 


5.3 


6.7 


4.8 


5.3 


16384 


2.4 


3.1 


2.0 


2.6 


32768 


1.0 


2.0 


1.5 


2.1 



0.1 



0.01 \r 

0.001 
1e-04 r 
1e-05 
1e-06 \r 
1e-07 



1 



^ 1 T-OAVAR 

— I F-OAVAR rough 

— I F-OAVAR without drift 

— I F-OAVAR circuiarized 



10 



100 1000 
integration time: i (s) 



10000 



100000 



Fig. 8. OAVAR for a White Frequency Noise sequence of length A'^ = 
65536 computed in the time domain and in the frequency domain, rough, 
after removing the linear frequency drift and after circularizing the sequence. 



F-OAVAR after drift removal /T-OAVAR 



the frequency domain variances we don't realize any extension 
of data manually as done in the computation of the Total 
variance. 

Table ( I VII ) compares the edf of the OAVAR for a Random 
Walk Frequency Noise computed after these processings. The 
best estimates are obtained by using the circularized sequence 
since the edf of the estimates are higher than for for the 
sequence after removing a linear frequency drift. Thus, the 
edf of the last estimate (r = N/2) is 2 times higher than the 
one of the estimate obtained in the time domain. This means 
that this estimate provides a noise level estimation V2 times 
sooner than the estimate computed in the time domain (e.g. 
265 days instead of 1 year) with the same accuracy. 

However, applying the circularization processing to another 
type of noise induced is a bias that has the same characteristic 
as a linear frequency drift on an Allan variance plot. Beside 
the behaviour characteristic of a white FM, figure [8] 
exhibits the r signature of a linear frequency drift in the 
Allan variance curve of the circularized sequence. Let us also 
notice the very long errorbars of the circularized sequence 
estimates. Therefore, the circularization process cannot be 
used in a real frequency deviation sequence which contains 
always different types of noise. Thus, we recommand to apply 
the spectral OAVAR over the residuals of a frequency deviation 
sequence, after removing the linear frequency drift. For a 
random walk FM, the estimate of OAVAR computed in the 
frequency domain after drift removal has an edf 1.5 times 
higher than the classical time domain OAVAR. It means that 
spectral OAVAR after drift removal is able to measure the 
random walk level of a sequence v^TTS times sooner than time 
OAVAR (e.g. 300 days instead of 1 year). 

Figure (|9]l compares the F-OAVAR variance computed after 
linear drift removal by least squares fit and the classical T- 
OAVAR variance. As shown in Table ( IVIb the upper plot shows 
that the edf of the F-OAVAR after drift removal for a Random 
Walk noise is less than the edf of the T-OAVAR for small r 
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— ■ — FLFM 
,9 1.5 - — « — RWFM 
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Fig. 9. Comparaison of the F-OAVAR computed after drift removal from 
noise sequences by least squares fit and the classical T-OAVAR for three noise 
types. The upper plot depicts the edf ratio and the lower plot depicts the bias. 
N = 56536. Results were obtained using Monte Carlo with 1000 trials. 

values. The lower plot shows that the F-OAVAR presents a bias 
of -10% for Random Walk noise. This bias can be explained 
by the fact the drift removal from a Random Walk sequance 
alters the spectrum of the noise at all the frequency values 
because a Random Walk contains a kind of linear drift feature 
intrinsicly. 

Let us remember that for a sequence without random walk 
FM (for atomic clocks), OAVAR computed in the frequency 
domain may be used directly and is more accurate than 
OAVAR computed in the time domain. 

B. OHVAR 

Figure ( fTOl l depicts the edf of the Overlapping Hadamard 
variance computed in the frequency domain F-OHVAR. It 
shows a very good agreement between the theoretical edf 
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F-OHVAR edf 



F-MAVAR vsT-MAVAR 




Averaging Factor 




Fig. 10. F-OHVAR edf for five noise types (a from -4 to 0) for sequences 
of length A'^ = 65536. The continous lines represent the theoretical edf 
computed by equation j96t . The symbols represent the edf obtained by Monte 
Carlo simulation with 1000 trials. 



Fig. 12. F-MAVAR and T-MAVAR for five noise types (a from -2 to +2) 
for sequences of length A'^ = 65536. The squares represent the F-MAVAR 
values and the dots represent the T-MAVAR values. Monte Carlo simulation 
with 1000 trials. 
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Fig. 11. Comparaison of the F-OHVAR and the T-OHVAR for five noise 
types. The upper plot depicts the edf ratio and the lower plot depicts the bias. 
A' = 56536. Results were obtained using Monte Carlo with 1000 trials. 

formula of equation ( |96] l and the edf obtained by Monte Carlo 
simulations. 

Figure (fTTl l shows that edf of the OHVAR estimator in 
the frequency domain is 2 to 4.5 higher than the edf of 
the classical OHVAR for the higher r = N/3 value. The 
lower plot depicts the bias defined by Bias = 100 x (1 — 
^F-OHVAR /T-OHVAR) . It is less than 10% for the five noise 
types and for all the r values. 

The Hadamard variance is not sensitive to linear frequency 
drifts. However, computing OHVAR in the frequency domain 
by using a FFT assumes also the periodicity of the sequence 
and may induce a high edge by connecting the last sample 
to the first one (see figure |7]- A). We performed then the same 
processings as previously in order to compare the effects of the 
drift removal and of the circularization of the sequence. For 



OHVAR also, the circularization should not be recommanded 
for processing frequency deviation sequences because it is only 
useful for noises with a < —2 and it degrades the variance 
estimates for the noises with a > —2 . On the other hand, the 
drift removal by substracting the best least squares line from 
the data gives good results for noises with a > — 2 . Hence, it 
is better to use the F-OHVAR directly without preprocessing 
in order to get better statistics than the T-OHVAR if the data 
does not contain a linear drift. 

C. MAVAR 

Figure (fT2l i shows a comparaison of the modified Allan 
variance computed in the frequency domain (F-MAVAR) and 
in the time domain (T-MAVAR) for five noise types with a 
from -2 to +2. We can notice clearly a huge bias of the F- 
MAVAR for a = +2 . 

For this reason, the use of MAVAR computed in the 
frequency domain should be avoided. 

VI. Conclusion 

We have presented a filter approach to analyze the different 
known frequency stability variances. Using this approach we 
derived formulae in the time domain identical to those known 
in the literature. We also demonstrated for the first time that the 
computation of these variances can be done in the frequency 
domain using a Discrete Fourier Transform of the studied 
signals. Such a computation provides estimates with better 
accuracy than the ones computed in the time domain, allowing 
the measurement of the low frequency noise levels sooner, 
i.e. with a shorter sequence. This advantage is particularly 
useful for studying the long term stability of atomic clocks. 
However, in the presence of linear drift, the periodicity of the 
sequence implicitely assumed by the use of the FFT algorithm 
may induce edges which degrade variance measurements if 
a random walk FM is present in the sequence. We have 
demonstrated that, in this case, we must first remove the linear 
frequency drift on a sequence before to compute a variance 
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in the frequency domain. Our work has proved that OAVAR 
computed in the frequency domain is the estimator which gives 
the quickest low frequency noise level (9 month instead of 1 
year). New estimators improving these characteristics with a 
more simple transfer function will be described in another 
paper [24]. 

Appendix : Equivalence of the Discrete-Time and 
THE Continuous-Time variances 

We have assumed T = 1 in (l44l i. Without this assumption 
the variance ct^ (m)could be written using dSTT) : 



(™) 



T 



T 



2r sm 



clvn? .Lj_ sin^ (tt/T) 



(98) 

Sy'{f)df. 



Using expression (|29] | of Sy{f) in ( l98T l we can write: 

■)2n r + Tr ( o;„2n+2 , 



{nfrnT) 



sin' (^/T) 



sin^ [^r 
[KT{f-nf,)f 



(99) 



df. 



The sine functions outside the sum sign are periodic, they 
can be passed inside the sum sign. Doing this and making the 
variable change i/ = / — nfs we can write: 



alim) 



22n +2^ f + TT-T 
" fe=-oo 

sin^ {-kTv) 



,2n+2 



{'KvmT) 



sin' {-nvT) 
dv (100) 



where we have interchanged the sum sign and the integration 
symbol. 

Equation ( II 00b simplifies to: 



22n .+oOgjj^2«+2 



^ -Sy {v)dv 



cr' (m) 

J-oo {■KmTv) 

= 4W(„)|r=mT. (101) 

The only difference between dlOlb and (l24l l is the integra- 
tion bounds. In equation (l24l) . S'y (/) is the single-sided PSD 
while Sy^ {u) in dlOll ) is the two-sided PSD. We conclude 
that ( |44] | and ( |24] | represent the same variance. 
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